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Abstract 

The angular momentum evolution of stars close to massive black holes (MBHs) is driven by secular torques. 
In contrast to two-body relaxation, where interactions between stars are incoherent, the resulting resonant 
relaxation (RR) process is characterized by coherence times of hundreds of orbital periods. In this paper, we 
show that all the statistical properties of RR can be reproduced in an autoregressive moving average (ARMA) 
model. We use the ARMA model, calibrated with extensive A-body simulations, to analyze the long-term 
evolution of stellar systems around MBHs with Monte Carlo simulations. 

We show that for a single-mass system in steady-state, a depression is carved out near an MBH as a result 
of tidal disruptions. Using Galactic center parameters, the extent of the depression is about 0.1 pc, of similar 
order to but less than the size of the observed "hole" in the distribution of bright late-type stars. We also find 
that the velocity vectors of stars around an MBH are locally not isotropic. In a second application, we evolve 
the highly eccentric orbits that result from the tidal disruption of binary stars, which are considered to be 
plausible precursors of the "S-stars" in the Galactic center. We find that RR predicts more highly eccentric 
(e > 0.9) S-star orbits than have been observed to date. 

Subject headings: black hole physics - celestial mechanics - Galaxy: center - stars: kinematics and dynamics 



1. introduction 

The gravitational potential near a massive black hole 
(MBH) is approximately equal to that of a Newtonian point 
particle. As a consequence, the orbits of stars are nearly Ke- 
plerian, and it is useful, both as a mental picture and as a 
computational device, to average the mass of the stars over 
their orbits and consider secular interactions between these el- 
lipses, rather than interactions between point particles. These 
stellar ellipses precess on timescales of many orbits, due to 
deviations from the Newtonian point particle approximation: 
there is typically an extended cluster of stars around the MBH, 
and there is precession due to general-relativistic (GR) ef- 
fects. Nevertheless, for timescales less than a precession time, 
torques between the orbi tal ellipses are coherent. 

It was first shown by iRauch & Tremaind (Q996) that such 
sustained coherent torques lead to much more rapid stochas- 
tic evolution of the angular momenta of the stars than normal 
relaxation dynamics. They called this process resonant relax- 
ation (RR). RR is potent ially important for a n umber of as- 
trophysical phenomena. Rauch & Ingalls] d 1998b showed that 
it increases the tidal disruption rate, although in their calcu- 
lations the effect was not very large since most tidally dis- 
rupted stars originated at distances that were too large for 
RR to be effective (in Section 17.21 we will revisit this ar- 
gument). By contrast with tidally disrupted main-sequence 
stars, inspiraling co mpact objects originate at dis tances much 
closer to the MBH dHopman & Alexanderil2005l) . The effect 
of RR on the rate of compact objects spiralling into MBHs to 
become gravitational wave so urces is therefore much larger 
( Hopman & Alexander 2006a). RR also plays a role in several 
proposed mechanisms for the origin of the "S-stars", a cluster 
of B -stars with randomized orbits i n the Galactic center (GC; 
e.g. iHopman & Alexander! I2006at iLevinl 120071: iPerets et alj 
20071 see Secfion l73l of this paper). 
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There have been several numerical studies of the RR 
process itself, which have verified the overall analytic 
predictions. H owever, since stellar orbits in A-body 
simul a tions (e.g. IRauch & Tremaind 1 1996b IRauch & Ingallsl 
| 1998t lAarsethl 120071: lHarfst et all 120081: fEilon et alj 120091: 
IPerets et al . 2009; Perets & Gualandris 2010) need to be in- 
tegrated for many precession times, the simulations are com- 
putationally demanding. All inquiries have thus far have been 
limited in integration time and/or particle number. In order 
to speed up the computation, several papers have made use of 
the picture described above, where the gravitational interac- 
tion bet ween massive wires are considered (Gauss's method , 
see e.g. IRauch & Tremaliiel 1 1996b iGtirkan & Hopmanl 120071: 
iTouma et al.ll2009h . 

It is common, when possible, to use the the Fokker- 
Planck formalism to carry out long-term simulations of 
the s tellar distribution in galactic nuclei dBahcall & Wol: 
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Hopman & Alex ander 2006a b). The current formalism how- 



ever is not directly applicable to the case when RR plays an 
important role. At the heart of all current Fokker-Plank ap- 
proaches is the assumption of a random-walk diffusion of an- 
gular momentum, whereas RR is a more complex relaxation 
mechanism based on persistent autocorrelations. 

In this work, we will show that the auto- regressive moving 
average (ARMA) model provides a faithful representation of 
all statistical properties of RR. This model, calibrated with 
special-purpose A-body simulations, allows us to carry out a 
study of the long-term effects of RR on the stellar cusp, thus 
far out of reach. 

The plan of the paper is as follows. In Section [2] we 
present an extensive suite of special purpose A-body sim- 
ulations, which exploit the near-Keplerian nature of stellar 
orbits and concentrate on the stochastic orbital evolution of 
several test stars. We use these simulations to statistically ex- 
amine the properties of RR for many secular timescales. In 
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Section [3] we introduce the ARMA model for the data anal- 
ysis of the A-body simulations. This description accommo- 
dates the random, non-secular (non-resonant) effects on very 
short timescales, the persistent autocorrelations for interme- 
diate timescales less than a precession time, and the random 
walk behavior for very long times (the RR regime). In Section 
2] we extend the ARMA model, using physical arguments, to 
a parameter space that is larger than that of the simulations. In 
Section [5] we calibrate the parameters using the results of the 
A-body simulations. Once the free parameters of the ARMA 
model are determined, we then use them in Section [6] as an 
input to Monte Carlo (MC) simulations which study the dis- 
tribution of stars near MBHs. In Section [7] we show that RR 
plays a major role in the global structure of the stellar dis- 
tribution near MBHs, and in particular for the young, massive 
B-type stars near the MBH (the "S-stars") in our GC. We sum- 
marize in Section[8] 

2. TV-BODY SIMULATIONS 

We have developed a special-purpose A-body code, de- 
signed to accurately integrate stellar orbits in near-Keplerian 
potentials; see Appendix [A] for a detailed description. We 
wrote this code specifically for the detailed study of RR. Such 
a study requires an integration scheme with an absence of spu- 
rious apsidal precession and one which is efficient enough to 
do many steps per stellar orbit while integrating many orbits 
to resolve a secular process. The main features of this code 
are the following: 

1. A separation between test particles and field particles. 
Field particles move on Kepler orbits, which precess 
due to their averaged potential and general relativity, 
and act as the mass that is responsible for dynamical 
evolution. Test particles are full A -body particles and 
serve as probes of this pot ential (Rauch & Tremaind 
[!9%tlRauch & Ingallslll998l) . 

2. The use of a four th-order mixed - variable symplectic 
(MVS ) integrator dYoshidalll990t I Wisdom & Holmanl 
I199U iKinoshita et al.1 1199 U ISaha & Tremaind 119921) . 
The MVS integrator switches between Cartesian coor- 
dinates (in which the perturbations to the orbit are cal- 
culated) to one based on Kepler elements (to calculate 
the Keplerian motion of the particle under the influence 
of central object). We make use of Kepler's equation 
(see, e.g., Danbv 1992) for the latter step. 

3. Adaptive time stepping and gravitational soften- 
ing. To resolve the periapsis of eccentric orbits 
(Rauch & Holman 1999) and close encounters between 
particles, we adapt the time steps of the p articles. We 
use the compact K2 kernel dDehne n 2001) for gravita- 
tional softening. 

4. Time symmetric algorithm. As MVS integrators gener- 
ally lose their symplectic properties if used with adap- 
tive time stepping, we enhance the energy conservation 
by time-symmetrizing the algorithm. 

Our code is efficient enough to study the evolution of energy 
and angular momentum of stars around MBHs for many pre- 
cession times, for a range of initial eccentricities. 

2.1. Model of Galactic Nucleus 



We base our Galactic nucleus model on a simplified GC 
templat^B It has three main components. (1) An MBH with 
mass M, = 4 x 10 6 M Q which remains at rest in the cen- 
ter of the coordinate system. (2) An embedded cluster of 
equal-mass field stars m = 10 M Q , distributed isotropically 
from 0.003 pc to 0.03 pc, which follow a power-law density 
profile n(r) oc r~ a . The outer r adius is chosen with refer- 
ence to lGiirkan & Hopmanl d2007), who show that stars with 
semi-major axes larger than the test star's apoapsis distance 
fapo = a(l + e) contribute a negligible amount of the net 
torque on the test star. The field stars move on precessing 
Kepler orbits, where the precession rate is determined by the 
smooth potential of the field stars themselves (see Appendix 
let and general relativity. The precession of the field stars is 
important to account for because for some eccentricities the 
precession rate of the test star is much lower than that of the 
"typical" field star, such that it is the precession of the latter 
that leads to decoherence of the system. The field stars do not 
interact with each other but they do interact with the test stars 
if the latter are assumed to be massive. In this way the field 
stars provide the potential of the cluster but are not used as 
dynamic tracers. (3) A number of test stars, used as probes 
of the background potential, that are either massless or have 
the same mass as the field stars, m = 10 Mq. We consider 
both massless and massive stars in order to st udy the effects 
of resonant friction (Rauch & Tremaine 1996). The test stars 
have semi-major axes of a — 0.01 pc and a specified initial 
eccentricity e. 

Following Equation (17) from Hopman & Alexander 
(2006 a|), who use the M m — a relation dFerrarese & M erritt 
2000; Gebhardt et al. 2000) which correlates the mass of a 
central black hole with the velocity dispersion of the host 
galaxy's bulge, we define the radius of influence as 

GM. ( M. V /2 

r " = ^ = 2 - 31pC UxlOSMj ' (1) 
The number of field stars within radius r is 

(\ 3— a 
— J , (2) 

where we assume that the mass in stars within equals that 
of the MBH, Nu = M ./m = 4 x 10 5 . We take a = 7/4 
(iBahcall & Wolifll976h . the classic result for the distribution 
of a single-mass population of stars around an MBH, which 
relies on the assumption that the mechanism through which 
stars exchange energy and angular momentum is dominated 
by two-body interactions. Hence the number of field stars 
within our model's outer radius is A(< 0.03 pc) = 1754. 
We summarize the potential for our Galactic nucleus model 
in TableQ] 

We evolve this model of a galactic nucleus for a 
wide range in eccentricity of the test stars, e = 
0.01, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 0.9 and 0.99. For each initial 
eccentricity, we follow the evolution of a total of 80 test stars, 
both massless and massive, in a galactic nucleus. Typically 
we use four realizations of the surrounding stellar cluster for 
each eccentricity, i.e., 20 test stars in each simulation. They 
have randomly oriented orbits with respect to one another, so 
that they experience different torques within the same clus- 
ter. The simulations are terminated after 6000 orbital periods 

1 Due to the effect of GR precession, the system is not scale-free and the 
masses need to be specified. 
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TABLE 1 
Galactic Nucleus Model 



Parameter 


Numerical Value 


M. 


4 X 10 b iW Q 


in 


10 Mq 


a 


1.75 


Th 


2.31 pc 


f*min 


0.003 pc 


?"max 


0.03 pc 


&t est star 


0.01 pc 


N(< r max ) 


1754 



(henceforth denoted by P = 2-ir^/a 3 /GM.) at a = 0.01 pc, 
roughly equivalent to three precession times for a star of ec- 
centricity e = 0.6, deep into the RR regime for all eccentrici- 
ties; see Section[5]for verification. Using our method of analy- 
sis, it would not be useful for our simulations to run for longer 
times as the stars would move significantly away from their 
initial eccentricities. In addition, the autocorrelation functions 
(ACFs) of their angular momentum changes drop to zero be- 
fore this time; see Figure [4] 

2.2. Illustrative simulations 

In Section [3] we present a description that captures all the 
relevant statistical properties of RR, and in particular has the 
correct autocorrelations. Here we consider several individ- 
ual simulations for illustration purposes, in order to highlight 
some interesting points and motivate our approach. We define 
energy E of the test star as 



E = 



GM. 

2a 1 



and dimensionless angular momenturrQ J as 



J = y/l - e 2 . 



(3) 



(4) 



Secular torques should affect the angular momentum evo- 
lution, but not the energy evolution of the stars (in the picture 
of interactions between massive ellipses of the introduction, 
the ellipses are fixed in space for times less than a preces- 
sion time, t < ip r cc, such that the potential and therefore 
the energy is time-independent). It is therefore of interest to 
compare the evolution of these two quantities. An example is 
shown in Figure Q] As expected, the angular momentum evo- 
lution is much faster than energy evolution; furthermore, the 
evolution of angular momentum is much less erratic, which 
visualizes the long coherences between the torques. 

In Figure|2]we show the eccentricity evolution for a sample 
population of stars with various initial orbital eccentricities. 
There is significant eccentricity evolution in most cases, but 
almost none for e = 0.99, the reasons for which we elucidate 
in Section[6] 

To quantify the rate at which the energy E and angular mo- 
mentum J change as a function of eccentricity e, we calcu- 
late the E and J relaxation timescales. We define these as 
the timescale for order of unity (circular angular momentum) 
changes in energy E (angular momentum J). We compute 

2 Throughout this paper, we use units in which angular momentum J and 
torque r are normalized by the circular angular momentum for a given semi- 
major axis, J c = \/ GM.a. All quantities are expressed per unit mass. 
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FIG. 1 . — Evolution of angular momentum J and energy E, normalized to 
the initial values, of a massless test star with semi-major axis a = 0.01 pc 
and eccentricity e = 0.6. Time is in units of the initial orbital period P (top 
axis shows time in Myr). The coherence of RR can be seen in the J evolution 
(precession time t prC c is ~ 2000 orbits), whereas the E evolution displays a 
much slower, more erratic, random walk. 



t |Myr| 




6000 



FIG. 2. — Eccentricity evolution of a sample of massless test stars in a 
fiducial galactic nucleus model. We show four realizations for five different 
initial eccentricities. Evolution is rapid for most eccentricities; however, the 
stars with the highest value e = 0.99 have sluggish eccentricity evolution. 



the following quantities: 



t E = 



E 2 



t.r 



<(A£) 2 

((aj) 2 ; 



-At, 



At, 



(5) 



(6) 



where AE and A J are the steps taken in a time At, and we 
take the mean ( ) over 80 test stars in each eccentricity bin. 

We plot the resulting timescales as a function of eccen- 
tricity in Figure [3] taking several At values to get order-of- 
magnitude estimates for these timescales. We find no signifi- 
cant difference between the results for the massless and mas- 
sive test stars. We note that t^ is only weakly dependent on e, 
whereas tj is of the same order of magnitude as t^ for e — > 
and e — > 1, but much smaller for intermediate eccentricities. 
It is in the latter regime that secular torques are dominating 
the evolution. In the following sections, we will give a de- 
tailed model for the evolution of angular momentum. Since 
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FIG. 3. — Timescale for order of unity (circular angular momentum) 
changes in energy E (angular momentum J) for massless test stars in our 
N-body simulations, presented on a log scale of time in units of Myr; see 
Equations j5J and (6). We find a weak trend for lower E relaxation times with 
increasing eccentricity (from ~ 400 Myr to 200 Myr). J relaxation occurs 
on very short timescales (~ 20 Myr) at eccentricities of 0.8 < e < 0.9. 
Stellar orbits with low eccentricities have very long J relaxational timescales 
(~ 1 Gyr) as the torque on a circular orbit drops to zero. This timescale 
increases again at very high eccentricities. 



the focus of this paper is mainly on secular dynamics, we do 
not further consider energy diffusion here, but in Appendix 
|B]we discuss the timescale for cusp formation due to energy 
evolution, and compare our results to other results in the liter- 
ature. 

3 . STATISTICAL DESCRIPTION OF RESONANT RELAXATION 

As a result of the autocorrelations in the changes of the an- 
gular momentum of a star, as exhibited in Figure Q] RR can- 
not be modeled as a random wal k for all times. This is also 
clear from physical arguments. Rauc h & Tremainel (119961) 
and several other papers take the approach of considering two 
regimes of evolution: for times t <C t picc , evolution is ap- 
proximately linear as the torques continue to point in the same 
direction. For times t t prcc , the torque autocorrelations 
vanish and it becomes possible to model the system as a ran- 
dom walk. Here we introduce a new approach, which unifies 
both regimes in a single description. This description is also 
useful as a way of quantifying the statistical properties of RR. 

We study the evolution of the angular momenta of the test 
stars using a time series of angular momenta at a regular spac- 
ing of one period. This choice is arbitrary, and below we show 
how our results generalize to any choice of time steps. The 
normalized ACF can be written as 



Pi 



((AJ s+t - (AJ))(AJ S - (AJ))) 
((AJ t -(AJ))2) 



(7) 



where A J s is the change in angular momentum between sub- 
sequent measurements at time s and it is assumed that the 
time series is stationary. A typical ACF from our simulations 
is shown in Figure [4] The fundamental feature of this ACF is 
that it is significantly larger than zero for hundreds of orbits, 
but with values much less than unity. 

Our goal in this section is to define a process that generates 
a time series that has the same ACF as that of the changes in 
angular momentum of a star. Such a time series has all the 
statistical properties that are relevant for the angular momen- 
tum evolution of the star, such as the same Fourier spectrum. 



1 1 Myr] 

0.0001 0.001 0.01 




10000 



FIG. 4. — Autocorrelation function for the differences A J in one particular 
simulation (noisy orange line) as a function of time in units of the initial 
orbital period P. The test star is massless and has initial eccentricity e = 0.6. 
The autocorrelation function is significantly larger than zero (if much less 
than one) for several hundred orbits. The smooth red line gives the theoretical 
ACF for an ARMA( 1 , 1) model with fa = 0.995, 6»i = -0.976. 



We have two motivations for doing this. First, defining such a 
process implies that once the parameters of the model are cal- 
ibrated by use of the A^-body simulations, RR is fully defined. 
Second, a generated time series with the same statistical prop- 
erties as the time series generated in the physical process can 
be used in MC simulations to solve for the long term evolution 
of RR. 

3.1. The autoregressive moving average model ARMA( 1,1) 

We now introduce the ARM A model, wh ich is often used 
in econometrics (see, e.g.. iHeij et al.l 120041) . We show that 
this model can generate time series that reproduce the ACF 
and variance of AJ t . The model is ad hoc in that it does 
not have a direct physical foundation. However, we will find 
physical interpretations for the free parameters of the model 
in Section S] We use a form of the model, ARMA(1,1), with 
one autoregressive parameter, <f>\, and one moving average 
parameter, 0\. The model can then be written as 



A x J t = (t>xA x J t -x 



9 ^ 



.(1) 



(8) 



The label "1" in this equation refers to the fact that the data 
have a regular spacing of one period. In this equation, cf>i and 
f?i are free parameters, and the random variable e is drawn 
from a normal distribution with 



< £ «> = 0; (eWeW) = a f §t 



(9) 



where o\ is a third free parameter and <5( S is the Rronecker 
delta. For such a model, the variance of the angular momen- 
tum step is 



(AJ, 2 ) = 



1 + 91 +26»ic 



1 



(10) 



and the ACF can be derived analytically (see Appendix iDl for 
details): 



Pt = <f> 



1 



i + (^i + 0i) 2 /(i-0i) 



(t > 0). 



(11) 
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From this expression, we see that the decay time of the ACF 
is determined by the parameter </>i, and B\ captures the mag- 
nitude of the ACF. As an example, we plot the ACF in Fig- 
ure|4]with parameters typical for RR. We emphasize the need 
for both an autoregressive parameter and a moving average 
parameter to reproduce the slow decay of the ACF. Neither 
parameters can replicate the features on their own. 

A useful reformulation of the ARMA(1, 1) model, which 
can be found using recursion, is 



AJ t = y^T/> fc e t -fc, 



(12) 



where 



4>k = (0i + 



h)<t>\- x 



0>0); ^o = L (13) 



Once the model parameters (<fii ,0i,a±) are found for a time 
step of one period, 5t = P, the variance after N steps can be 
computed. The displacement after N steps is given by 



N 



N 



(14) 



n=l k=0 



where (n — k) > 0. We denote the variance of the displace- 
ment after N steps by 



(15) 



such that 



N oo N 



n=l k=0 m=l 1=0 
N oo N oo 

= CT 1 E E E E V'fcV^ n—k,m—l i (16) 
n=l k=0 m—l 1=0 

which can be decomposed as 

JV m-l AT-1 

m — n 

+ MN -k) 

m—2 n— 1 k—1 
oo N min(/c+m— l.AT) 



k—1 m—l 



(17) 



This expression can be cast in a form without summations by 
repeatedly using the properties of the independent normal ran- 
dom variables e t in Equation (O, resulting after some algebra 
in 



V {N)o^=N+^±p- 

1-01 



+ 



2AT-1 + 

2 

'jV-2 



20f - 0! - 1 

1-01 
0i(i-O n 



1 



(18) 

Note that in the special case that N = 1, Equation ( flOl l is 
recovered. We plot this expression in Figure [5] 

Three characteristic timescales are clearly discernible in 
this figure. On a short timescale, the angular momentum 
diffusion is dominated by the usual two-body gravitational 
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FIG . 5. — Variance in angular momentum J after time t (V(t); from Equa- 
tion )18t ) as a function of time in units of the initial orbital period P, with 
(pi = 0.999, 0i = -0.988 and a\ = 10 -4 . Note that there is first a linear 
part ("NR"), then a quadratic part ("coherent torque"), and then again a linear 
part ("RR"). 



scattering (i.e., non-resonant relaxation, hereafter NR). In this 
regime the angular momentum deviation scales as ^/t. On an 
intermediate timescale the orbital angular momentum evolves 
due to the nearly constant orbit-averaged secular torques from 
the other stellar orbits, and thus the angular momentum de- 
viation scales linearly with time. On a yet longer timescale 
the secular torques fluctuate due to orbital precession, and the 
V< scaling is reinstated, albeit with a much higher effective 
diffusion coefficient. Thus the ARMA-driven evolution is in 
agreement with what is expect ed for a combinatio n of RR and 
NR acting togeth er; see also iRauch & T remaine (199g) and 
Eilon et al.l (120091) for discussion. 

For the typical values of <\>\ and Q\ that we find, these differ- 
ent regimes can be identified in the expression for the variance 
V(N); Equation (1 181 . For very short times, the expression is 
dominated by the first term, which can thus be associated with 
NR. The second term and third term are both proportional to 
A^ 2 for short times (this can be verified by Taylor expansions 
assuming N(pi 1; the linear terms cancel), and as such rep- 
resent coherent torques; and they become proportional to N 
for large time, representing RR. In Figure[5]we plot Equation 
( fT8l for typical values of the parameters of the model. 

Equation ( [18] ) is useful for finding the variance after N pe- 
riods, but most importantly to define a new process with ar- 
bitrary time steps. We can now define a new ARMA(1, 1) 
process with larger time steps, but with the same statistical 
properties. If the new process makes steps of N orbits each, 
we write the new model as 



A N J t = (f> N A N J t ^ N + 6>Are|^ + e\ N '; 



(19) 



<[ e (A °] 2 > = a^, 
where the parameters [cj)(5t),9(St)] are related to (</>i, 8\) b)0 

4>{8t) = 4 t/P ; 9(St) = -[-8^/?, (20) 
and in particular, the ARMA parameters for steps of N orbits 

3 For cf>i we use the fact that the theoretical ACF of an ARMA(1, 1) model 
decays on a timescale P/ In rf>± . We have no interpretation for the functional 
form of 9\ , which was found after some experimenting. We confirm nu- 
merically that this form leads to excellent agreement with simulations with 
arbitrary time; see Figure|6] 



6 



MADIGAN, HOPMAN AND LEVIN 




FIG. 6. — Monte Carlo simulation of RR with the model parameters 
tjn = 0.999, 6»i = -0.988, and ai = 10~ 6 (chosen for purely illustra- 
tive purposes), and the time steps in units of the orbital period indicated in 
the legend. After 10 4 orbits, the distribution is very similar for all time steps 
chosen, confirming that the diffusion in J is independent of the time step. 



<t> N = 4>{NP) = 4> 



= 0(NP) = -[-6»i] Jv . (21) 



Since after N small steps of one period, the variance should 
be the same as after 1 single step of time N orbits, we can find 
the new a to be (see Equations ( [Tol l. ( TToT l) 



1 



>N 



1 + 6% + 20 



nun 



-V(N). 



(22) 



Once the ARMA parameters for steps of one orbital period 
are found, they can be generalized to a new ARMA model via 
Equation ( fT9l ), so that the time step equals an arbitrary num- 
ber of periods (not necessarily an integer number), and the 
new parameters are found by the use of Equations (OH and 
( f22b . We will exploit these relations in MC simulations (Sec- 
tion |6j, in which we use adaptive time steps. We have tested 
the equivalence of the models by considering an ARMA dif- 
fusion process over a fixed time interval but with varied time 
steps. We show an example in Figure [6] where we have taken 
a delta function for the initial ./-distribution, and have plot- 
ted the resulting final distributions after 10 4 orbits for sev- 
eral different time steps which vary by two orders of mag- 
nitude. We find excellent agreement between these distribu- 
tion functions which verify the generalization of (<f>i , Q\ , o\ ) 
to (4> N ,9 N ,a N ). 

4. INTERPRETATION OF THE ARMA PARAMETERS AND 
EXTENSION OF PARAMETER SPACE 

The ARMA(1, 1) model presented here contains the three 
free parameters (0i, 9%, ax). In this section, we interpret these 
parameters in the context of a stylized model of a galactic nu- 
cleus. Our purpose in doing so is twofold. First, it relates the 
ARMA model to the relevant physical processes that play a 
part in the angular momentum evolution of stars near MBHs, 
and gives insight into how the parameters depend on the stel- 
lar orbits, in particular on their eccentricity. This allows us to 
define, and discuss, RR in the language of the ARMA model. 
Second, in Section [5] we calibrate the parameters as a func- 
tion of eccentricity using iV-body simulati ons for one spe- 
cific configuration — as specified in Section l2~T1 In order for 
the parameters to be also valid for other conditions than those 



studied, one needs to understand how they vary as a function 
of the quantities that define the system. The model presented 
in this section, which is based on physical arguments, allows 
us to generalize our results to models we have not simulated 
(models with different stellar masses, semi-major axes, etc.) 
The arguments presented here do not aspire to give a full the- 
oretical model of RR, but rather parameterize the model in a 
simplified way that is useful for generalization to other sys- 
tems. We will use this model in our MC simulations (Section 

4.1. Non-resonant relaxation (NR): The parameter a 

For very short times NR dominates the angular momentum 
evolution of a star and the variance of J is expected to be of 
order -P/£nr after one period, where 



^nr — A 



NR 



i_ 

m ) N < In A 



P 



(23) 



is the NR time dRauch & Tremainel fl996i). defined as the 
timescale over which a star changes its angular momentum 
by an order of the circular angular momentum at that radius. 
In this expression, N < is the number of stars within the ra- 
dius equal to the semi-major axis of the star, A = M./m, 
and Anr = 0.26, chosen to match the value of ts from the 
iV-body simulations; see Figure[3] Since this variance should 
be approximately equal to a\, it follows that 



! ° 1 M. 



A^lnA 



.4 



(24) 



NR 



Here f a is a (new) free parameter which, calibrated against 
our iV-body simulations, we expect to be of order unity (the 
deviation from unity will quantify a difference in and NR 
component of tj). In our iV-body simulations we find that o\ 
is an increasing function of e; see Figures [12] and [T4j This 
suggests that NR is also an increasing function of e which we 
attribute to the increase in stellar density at small radii which 
a star on a highly eccentric orbit passes through at periapsis. 
To account for this in our theory we fit f a as a function of e. 

4.2. Persistence of coherent torques: The parameter (j) 

We now define a new timescale over which interactions be- 
tween stars remain coherent, the coherence time t^,, which 
corresponds to the time over which secular torques between 
orbits remain constant. This is the minimum between the pre- 
cession time of the test star and the median precession time 
of the field stars. The latter timescale is of importance as 
there exists at every radius an orbital eccentricity which has 
equal, but oppositely directed, Newtonian precession due to 
the potential of the stellar cluster and GR precession. A star 
with this particular orbital eccentricity remains almost fixed 
in space; for our chosen model of a galactic nucleus the ec- 
centricity at which this occurs is e ~ 0.92 at 0.01 pc. Its 
coherence time, however, does not tend to infinity, as the sur- 
rounding stellar orbits within the cluster continue to spatially 
randomize. 

The timescale for the Newtonian precession (i.e., the time 
it takes for the orbital periapsis to precess by 2ir radians) is 
given by 



tt(2 - a) 



N^rr, 



P(a)f(e,a)- 



(25) 
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(see Appendix ICl for derivation), where A< is the number of 
stars within the semi-major axis a of the test star, P(a) its 
period, and 



/(e,a = 7/4)- 1 
= 0.975(1 -e)- 1/2 + 0.362(1 - e) 



0.689 



(26) 



The timescale for precession due to general relativity is given 
by 



j-GR _ 1 M „2\ ac 



t = 2(1 - e 2 V 

''prcc o V / , 



-P(a) 



proc g v- - 

(Einstein 1916). The combined precession time is 

1 1 



(27) 



(28) 



As stated, the coherence time for a particular star is the min- 
imum of its precession time and that of the median precession 
time of the surrounding stars, 



t<j> = U min [iprecK e), t prec (a, e)] 



(29) 



where e is the median value of the eccentricity of the field 
stars (so e = y/1/2 for a thermal distribution) and is a 
second (new) free parameter. 

We now have all the ingredients to define the RR time. This 
is the timescale over which the angular momentum of a star 
changes by order of the circular momentum. The angular mo- 
mentum evolution is coherent over a time of order t$, and the 
change in angular momentum during that time is ri^, where r 
is the torque exerted on the star's orbit. N ormalized to the cir- 
cular angular momentum, this torque is dRauch & Tremainq 
fl996h 



M. P 



(30) 



w here the linear dependence on eccentricity was determined 
bv lGurkan & Hopmanl (12007). We adopt their value of A T — 
1.57 which was determined for a different stellar mass profile 
m = IMq, a = 1.4 where M. = 3 x 10 6 M Q ,r h = 2pc. 
Assuming that the evolution is random for longer times, this 
leads to the definition 



iRR(.E, J)=[^j- 
Tt<k 



(31) 



The theoretical ACF of an ARMA(1, 1) model in Equation 
( [TO shows that it decays on a timescale Pj \n4>\. Hence we 
find that scales as exp(— P/t$). However, even if the co- 
herence time is very long, the evolution of the orbit may not 
be driven by secular dynamics. The secular torques on an or- 
bit are proportional to its eccentricity so if the eccentricity is 
very small secular effects are weak. As a result, the evolution 
is dominated by two-body interactions, which have a vanish- 
ingly short coherence time. Within our formalism, this effect 
can be accounted for by multiplying t$ by a function S, 



exp 



P 

Std 



(32) 
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FIG. 7. — Critical eccentricity e c[lt for which NR and RR mechanisms 
work on comparable timescales, as a function of semi-major axis a from the 
MBH in our Galactic nucleus model. The shaded region in this plot indicates 
the eccentricities for which RR is the dominant mechanism. The lower values 
of e cr it decrease with distance to the MBH as RR becomes more effective (r 
and tj, increase with decreasing a). At a ~ 0.004 pc this value rises again 

in our model as t^R, decreases rapidly as r" 1 / 4 . The higher values of e cr it 
are due to GR effects which decrease tj, toward the MBH and substantially 
reduce the effectiveness of RR. The lower values of e cr it are due to the fact 
that t oc e. 



which interpolates between the secular and two-body regime. 
We remind the reader that the subscript "1" is used to em- 
phasize that this is the value of cf> for a time step of one or- 
bital period. The exponential transition between NR and RR 
regimes is somewhat arbitrarily chosen but fits the data rea- 
sonably well, k > is a new parameter which determines the 
steepness of the transition (a third free parameter), and e cr it 
is defined as the critical eccentricity at which the NR and RR 
times are equal (Equations ( 1231 , (TJlTl): 



t(a,e) 



In A 



1/2 



(34) 



with 



Inserting relevant values into the above equation and solv- 
ing numerically using a Newton-Raphson algorithm, we cal- 
culate e cr it ~ 0.3 at a — 0.01 pc which is in good agreement 
with the A^-body simulations. We plot e cr ; t as a function of a 
in Figure|7] RR becomes more effective at lower a as the mass 
precession time increases, and the value of e cr it decreases. At 
the lowest a there are two roots to Equation d34b : this is due to 
GR effects which compete with mass precession and decrease 
the coherence time to such an extent as to nullify the effects 
of RR at high eccentricities. These results are derived for a 
simplified galactic nucleus model and will change for differ- 
ent assumptions on the mass distribution; nevertheless we find 
the concept of e cr it useful for understanding the relative im- 
portance of RR and NR at different radii. 

4.3. Magnitude of resonant relaxation steps: The parameter 

9 

For times less than t<j>, there are coherent torques between 
stellar orbits; see Equation d30T >. The expected variance after 
a time t<j> is then 
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Alternatively, we can use the ACF (see Appendix [D] for 
derivation) to find 



(AJ|) = of 



Equating Equations d35l > and d36l > gives two values for 6\ 



(36) 




is +n 



4(1 - 2 )r 2 P 2 



(37) 




exp 



where fg is a fourth free parameter. 



4(1 - 2 )r 2 P 2 



Taking the positive root we find an excellent match to the data. 
For values of e > e cr ;t the first term + 4>i) ~ 1 ar, d 
recognizing the Taylor expansion of — exp(— x) « x — 1 we 
simplify this expression to 



(38) 



We summarize this section by clarifying that we have four 
free parameters in our ARMA model (fg, f$, f a , k) fully de- 
terminable by our TV-body simulations, with which we cali- 
brate our model for use in MC simulations. 

5. RESULTS: ARMA ANALYSIS OF THE TV-BODY SIMULATIONS 

We now return to the TV-body simulations, and use the time 
series of angular momenta that are generated to calibrate the 
ARMA model. These parameters fully define both the RR 
and NR processes. Once they are known, they can be used to 
generate new angular momentum series that have the correct 
statistical properties, much like what is done in regular MC 
simulations. We will exploit this method in Section|6] 

When a test star has mass, there will be a back-reaction on 
the field stars kn own as resonant friction, analogous to dy- 
namical friction dRauch & Tremainel[T99 6). In our simula- 
tions, we consider both the case that the test star is massless 
and that it has the same mass as the field stars to calibrate the 
model parameters. Differences in the parameters then result 
from resonant friction. 

We use the "R" language for statistical computing to 
calculate the ARMA model parameters from our simulations. 
We input the individual angular momentum time series of a 
test star, and use a maximum likelihood method to calculate 
(pi, d\ and u\. We find that the drift term (A J) for stars of all 
eccentricities is 0(1O~ 6 ). 

In Figure [8] and [9] we show scatter plots of the model pa- 
rameters 1 — 01 and for simulations in which the test 
stars are massless and massive respectively. We present these 
quantities rather than <j>\ and 9\ themselves because they span 
several orders of magnitude and, physically speaking, the rel- 
evant question is by how much the parameters differ from (mi- 
nus) unity. For intermediate eccentricities (0.4 < e < 0.9), 
1 — 0i is much smaller than unity, indicating that the an- 
gular momentum autocorrelations are indeed persistent, and 
there are significant torques between stellar orbits. However, 
the coherent tor ques are mixed with NR no ise (the classical 
NR as treated in iBinnev & Tremaind (120081) ). and as a result 
1 + 0i <C 1 as well. These values for 1 — <f>\ and 1 + 6\ con- 
firm that although the ACF for RR is finite for long times, it is 



e=0.01 
e=0.1 
e=0.2 
e=0.3 
e=0.4 
e=0.6 
e=0.8 
e=0.9 
e=0.99 




0.001 



0.01 



0.1 



FIG. 8. — Scatter plot of (1 — </>i) vs (1 + £>i) as found from data-analysis 
on the TV-body runs, for the case of massless test stars. For most cases, the 
parameters cluster, leading to a concentration of points around small values; 
reasons for exceptions are discussed in the text. From inspection of several 
individual cases, we see that the test stars of intermediate eccentricities which 
have values tj>i, \B\ | < 1 (to the bottom left of plot) have non-zero autocorre- 
lations for long times, similar to Figure|4] This is as expected from Equation 
{TT}. Test stars with <f>i,6i ~ have ACFs that are close to zero everywhere, 
even for short times. 
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FIG. 9. — Scatter plot of (1 — (pi) vs (1 + 0\) as found from data analysis 
on the TV-body runs, for the case of massive (m = IOMq) test stars. We 
find that, for most eccentricities, the values for 1 — <pi and 1 + 9\ cluster 
around those for the massless cases. 



much smaller than unity — which would be the case if there 
was no NR noise at all. 

For very small and very large eccentricities, we find 0i w 
6*i w 0, so there are no persistent torques of significance. The 
reason for the absence of coherent torques on the high and low 
eccentricity limits are different: for very small eccentricities, 
the secular torque on a ste llar orbit approaches zero as t oc e 
dGiirkan & Hopmanl2007h . At the high eccentricity end, there 
are significant torques, but the coherence time is so short due 
to general relativity, that the persistence of these torques is 
neglig ible, and there is no c oherent effect (quenching of RR: 
see lMerritt & Vasilievll20 1 Oh . Test stars with eccentricities of 
0.2 < e < 0.4 have a large scatter in their 1 — 01 and 1 + &i 
values. This is due to their proximity to e cr it, the eccentricity 
at which NR and RR compete for dominance at this radius in 
our model. 

In Figures [10] and [TT] we compare the median values of 
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FIG . 10. — Median value of 1 — 0i , plotted on a log scale, for both massless 
and massive runs, computed for 80 test stars in each eccentricity bin. Closed 
boxes denote the 45th and 55th percentiles, while the lines indicate the one 
sigma values. Note that for very high eccentricities the values increase, which 
in the interpretation of Section|4]means that the coherence time has decreased. 
This is consistent with the fact that for e > 0.92, GR precession becomes 
dominant over mass precession. The deviation at low eccentricities between 
values for massless and massive particles may due to resonant friction. 
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FIG. 1 1. — Median value of 1 + 0\, plotted on a log scale, for massless 
and massive runs, computed for 80 test stars in each eccentricity bin. Closed 
boxes denote the 45th and 55th percentiles, while the lines indicate the one 
sigma values. The results for massless and massive particles follow each 
other closely except at the low-eccentricity end where there is a large scatter 
in the data. 



1 — 01 and l + for both massive and massless test stars. 
There is significant scatter in the data near 0.2 < e < 0.4, 
where NR relaxation competes with RR, as can be seen from 
the 45th, 55th and one sigma (34.1 — 84.1) percentiles. The 
differences between the median values for massless and mas- 
sive test stars are very small, such that we did not see a strong 
feedback effect due to resonant friction, except for small stel- 
lar eccentricities. Resonant friction appears to drag stars away 
from circular angular momentum. The same deviation is seen 
for the a i model parameter in Figure [121 

From here on, we use only the results of the massless simu- 
lations as we have a larger range of test star eccentricities and 
have found the ARMA parameters {<f>\, 6%, a%) to be similar 
in value. We refer the reader to Table|2]for exact quantities. 

The ARMA model parameters for both very low eccentric- 
ities (e ~ 0.01) and high eccentricities (e > 0.4) tend to 



FIG. 12. — Data for the model parameter a\ as a function of eccentricity for 
massless and massive test stars. Plotted are the 45th to 55th percentiles (open 
boxes), one sigma values (lines), and 50th percentile for each eccentricity. 
The results for massive test stars are shifted along the a>axis for clarity. The 
values follow each other closely with the exception of low eccentricities. 



TABLE 2 

Median values of 1 — fa, 1 + 6>i and ai for massless test stars 

AS A FUNCTION OF ECCENTRICITY. 



e 


I- fa 


l + 6»i 


erj (10- 4 ) 


0.01 


0.992 


1.004 


1.853 


0.1 


0.993 


1.002 


1.828 


0.2 


0.991 


1.003 


1.866 


0.3 


0.569 


0.609 


1.923 


0.4 


0.006 


0.024 


2.027 


0.6 


0.005 


0.025 


2.184 


0.8 


0.004 


0.031 


2.379 


0.9 


0.004 


0.033 


2.460 


0.99 


0.047 


0.156 


2.626 



cluster together. For these cases one form of relaxation, ei- 
ther NR or RR, dominates the form of the ACFs, and taking 
mean values of (<^>i,#i) gives an accurate representation of 
the population. However, near the e cr i t boundary where the 
two relaxation effects compete for dominance, (<^i, 9\) do not 
consistently cluster but rather choose NR or RR values which 
can vary significantly. Hence we choose the median value in 
each case to compare with our theory. 

In Figures[T3]and[T4]we compare the experimental median 
values of (1 — <p\), (1 + 0\) and cri to our theoretical model. 
We find that good agreement between the iV-body simula- 
tions and the model is obtained by assuming values for the 
free parameters (f^, fg, k) = (0.105, 1.2, 30). The value of 
fcf, shows that the coherence time t$ is much less than the pre- 
cession time t r 



of e we find 
U = 0-52 



prec , see Equation ( |29l . Fitting f a as a function 
0.62e - 0.36e 2 + 0.21e 3 - 0.29^- (39) 



6. MONTE CARLO SIMULATIONS AND APPLICATIONS 

Our aims in this section are to (1) explore the long-term 
evolution of stars around an MBH, (2) investigate the deple- 
tion of stars around an MBH due to RR in context of the 
"hole" in late-type stars in the GC, and (3) study the evolu- 
tion of the orbits of possible S-star progenitors from different 
initial setups. Direct iV-body simulations of secular dynamics 
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FIG. 13. — Media n value of (1 — <j>i) and (l + 6>i), compared with the theo- 
retical values from B2II381 . There is a sharp transition between eccentricities 
0-3 < e < 0.4 where the two effects of NR and RR compete. A second 
transition occurs at high eccentricities (e > 0.92) as the coherence time, and 
hence <f>i, decreases due to general relativistic effects. The theory does not 
exactly match the median value of the data at 0.2 < e < 0.3. However we 
can vary the value of e cr ; t to reconcile this difference and we find that this is 
not important for our results in the next section. 
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FIG. 14. — Median values of the model parameter ai (related to non- 
resonant relaxation) and theoretical fit, as a function of eccentricity for mass- 
less test stars. 



are challenging due to the inherent large range in timescales. 
The precession time is orders of magnitude larger than the or- 
bital time, and in order to study relaxation to a steady-state, 
the system needs to evolve for a large number of precession 
times. Currently it is not feasible to do this using A-body 
techniques. Instead we draw on MC simulations to study 
long term secular evolution. For this purpose, we have de- 
veloped a two-dimensional (energy and angular momentum) 
MC code. The main new aspect compared to earlier work is 
that it implicitly includes evolution due to RR. The angular 
momentum diffusion is based on the ARMA(1, 1) model pre- 
sented in Section|3] For the en ergy evolution we use a similar 
method as in IHopmanl (120091) . which was in turn based on 
Shapiro & Marchant (1978). A key feature of this method is 
a cloning scheme, which allows one to resolve the distribu- 
tion function over many orders of magnitude in energy, even 
though the number of particles per unit energy is typically a 
steep power-law, dN/dE oc _E -9 / 4 . We will not describe 
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FIG. 15. — RR and NR timescales in units of Myr as a function of semi- 
major axis a in our Galactic nucleus model for different eccentricities. 

the cloning scheme here; f or d etails we refer the reader to 
Shapiro & Mar chantl (119781) and IHopmanl (120091) . 

We base our model of a galactic nucleus on parameters rel- 
evant for the GC; see Section |2] We set the MBH mass as 
M, = 4 x 10 6 M©, which affects the physics through the GR 
precession rate and the loss cone. We linearize the system by 
considering the evolution of test stars on a fixed background 
stellar potential. We assume a single-mass distribution of stars 
of either m = IMq or m = IOMq, with the number of stars 
within a radius r, N{< r), following a power-law profile as 
in Equation (|2). 

In Figure Q3] we plot the NR and RR timescales for stars 
with different orbital eccentricities as a function of semi- 
major axis for this model with m = IOMq. A cuspy min- 
imum arises in the RR timescales where GR precession can- 
cels with Newtonian mas s precessio n. In their GC m odel 
iKocsis & Tremainel (120101) find, using lEilon et ail (120091) pa- 
rameters, the cuspy minimum of RR near r = 0.007 pc (see 
their Figure 1), which falls precisely within the range of val- 
ues in our model (0.004 — 0.02 pc). Figure [T6l shows the RR 
time as a function of eccentricity e for different stellar masses 
within 0.01 pc. Mode ls with different enclosed masses will 
be used in Section lTTSI 

In reality, several stellar populations are present in the GC, 
and our choice for a single mass model is therefore a sim- 
plification. The slope a = 1.75 is smaller than would be 
expected for strong mass -segregation (lAlexander & Hopmanl 
2009t iKeshet etlil 2009). but was chosen to make the model 
self-consistent as a collection of intera cting single mass parti - 
cles will relax to have this distribution dBahcall & Wolfll976l) . 

6.1. Angular momentum 

For the angular momentum evolution of a test star, we go 
through the following steps: we first calculate the ARM A 
model parameters for one orbital period (cj)-y,6\, cr\) through 
Equations (d32l. d38l . d24"ii). We then find the model parame- 
ters (4>n, Qn, ctjv) for the time step St, where N(= St/ P) e 
R > 0, using Equations (f2Qb - d22l . Finally, the step in angu- 
lar momentum is given by Equation ( fl9] > 



K A T J-f) *W J-^W. 



(40) 



•N- 



6.2. Energy 
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FIG. 16. — Resonant relaxation timescale t rr for varying stellar mass (Nm 
where -m = IOMq) within 0.01 pc as afunction of e. In this plot the density 
power-law index a = 1.75. The e at which RR is most efficient decreases 
with decreasing mass. The solid black line indicates the value of t rr if the 
effects of GR precession are not taken into account (and is independent of 
number N). In our fiducial model the optimal efficiency (shortest t Ir ) is for 
Nm = IOOOMq. At this radius, GR effects are important for RR and t rr is 
not independent of N. 

The energy of a star is given in units of the square of the 
velocity dispersion at the radius of influence, vl, such that 
E = r%j1a. The step in energy during a time St is 



AE = £E 



St \ 



1/2 



(41) 



where £ is an independent normal random variable with zero 
mean and unit variance, and £nr is given by Equation d23l 
with ^ NR = 0.26. 

6.3. Initial conditions and boundary conditions 

For our steady-state simulations, we initialize stars with 
energies E = 1 and angular momenta drawn from a ther- 
mal distribution, dN(J)/dJ oc J. We follow the dynamics 
of stars between energy boundaries 0.5 = E m i n < E < 
£ max = 10 4 (for our model this corresponds to 2.31 pc > 
a > 1.155 x 10~ 4 pc). Stars that diffuse to E < E min are 
reinitia lized, but their time i s not s et to zero (zero-flow so- 
lution; |BaEcalL&\^[l97l [l977|). Stars with E > E max 
are disrupted by the MBH and also reinitialized while keep- 
ing their time. Stars have angular momenta in the range 
Jic < J < 1, where J\ c is the angular momentum of the last 
stable orbit. Stars with J < J\ c are disrupted and reinitialized, 
while stars with J > 1 are reflected such that J — >• 2 — J. 

6.4. Time steps 

In the coherent regime, the change in angular mo- 
mentum durin g a time St is of order SJ ~ rSt ~ 
A T e(m/M,)\/ N < (St / P) . The changes in angular momen- 
tum should be small compared to either the size of the loss 
cone Jic, or the separation between J and the loss cone, 
J — Jj c ; and compared to the distance between the circular 
angular momentum J c and J, i.e., compared to 1 — J. We 
therefore define the time step to be 

St = ft r(o, e)~ x min [J c - J, max (| J - J\ c \,Ji c )] , (42) 

where / 4 < 1. After some experimenting we found that the 
simulations converge for f t < 10~ 3 , which is the value we 
then used in our simulations. 



6.5. Treatment of the loss cone 

Stars with mass m and radius i?* on orbits which pass 
through the tidal radius 



n 



( 2M. 

\ m 



1/3 



(43) 



are disrupted by the MBH. The loss cone is the region in an- 
gular momentum space which delineates these orbits and is 
given by 

Ji c = y/2GM.r t . (44) 

Our pr escription for the time step , which is similar to that 
used in lShapiro & Marchantl (119781) . ensures that stars always 
diffuse into the loss cone and cannot jump over it. Stars are 
disrupted by the MBH when their orbital angular momentum 
is smaller than the loss cone, and they pass through periap- 
sis. Only if floor(i/P) - floor [{t - St) / P] > 0, do we con- 
sider the star to have crossed periapsis during the last time 
step and hence to be tidally disrupted (or directly consumed 
by the MBH in case of a compact remnant). In that case we 
record the energy at which the star was disrupted and initial- 
ize a new star, which starts at the time t at which the star was 
disrupted. 

7. RESULTS 

7.1. Evolution to steady-state 

We first consider a theoretical benchmark problem, that of 
the steady-state distribution function of a single-mass popula- 
tion of stars around an MBH. With m = 1M Q this solution 
appears after ~ 100 Gyr (10 energy diffusion timescales; see 
Appendix IB1. In order to highlight the differences caused by 
RR relative to NR, we define the function 



9{e,j 2 ) 



E 5 ' A d 2 N(E,J 2 ) 



(45) 



J 2 d\nEd\nJ 2 ' 

For a iBahcall & Woil (119761) distribution without RR, this 
function is constant. In lToonen et alJ (120091) . simulations sim- 
ilar to those presented here were used, except that angular 
momentum relaxation was assumed to be NR. It was shown 
that indeed g(E, J 2 ) is approximately constant for all (E, J) 
for the NR case; see Figure Al in that paper. 

In Figure [T7] we show a density plot of the function 
g(E, J 2 ) from our steady-state simulations, which illustrates 
several interesting effects of RR. At low energies, g(E, J 2 ) 
is constant which is to be expected since the precession time 
is of similar magnitude as the orbital period, such that there 
are no coherent torques. At higher e nergies (E > 10), there 
are far fewer stars than for a classical Bahca ll & Wolii (119761) 
cusp. This is a consequence of enhanced angular momen- 
tum relaxation, where stars are driven to the loss cone at a 
higher rate than can be replenished by energy diffusion, and 
was a nticipated by one-dimension al Fokker-Planck calcula- 
tions dHopman & Alexand er 2006a, see their Figure 2). At the 
highest energies considered (E > 100) most orbits accumu- 
late close to the loss cone, and there are a few orbits populated 
at larger angular momenta, with a "desert" in between. 

To further illustrate the reaction of the stellar orbits to RR, 
in Figure [18] we show the distribution of eccentricities of stars 
in slices of energy space. As in Figure [17] we see that for 
higher energies the distribution is double peaked, with most 
stars accumulating at the highest eccentricities, and another 
peak at eccentricities around e = 0.2. This distribution can 
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FIG. 17. — Two-di mensional distrib ution of stars, normalized such that for 
an isothermal Bahca ll & Woll 1 1976) profile the distribution should be con- 
stant (see Equation |45j). Note that angular momentum is in units of the 
circular angular momentum, so that the loss cone (indicated by the red line) 
is not constant. At high energies it is clear that there is a depletion of stars, 
and that the distribution of angular momenta is far from isotropic. 




FIG. 18. — Eccentricity distribution at fixed energies. For high energies the 
distribution of eccentricities is bimodal due to RR, whereas for low energies 
the distribution is isotropic with a cutoff at the loss cone. 



be understood as follows: RR is very effective at intermedi- 
ate eccentricities, where torques are strong and the coherence 
time is very long. As a result, the evolution of such eccentric- 
ities occurs on a short timescale. At high eccentricities, the 
coherence time is short due to GR precession, while at low 
eccentricities the torques are very weak. At such eccentrici- 
ties, evolution is very slow. The eccentricity therefore tends 
to "stick" at those values, leading to the two peaks in the dis- 
tribution. It is in teresting to note that A-body simulations by 
iRauch & Ingallsl ( 1199 8) also show a rise of the angular mo- 
mentum distribution near the loss cone. This result however 
may be affe cted by the new d ynamical mechanism recently 
described bv lMerritt et al.l (1201 ll) . 

In Figure[l9]we show the rate of direct capt ures of stars by 
the MBH per log energy. Loss cone theory (|Fr ank & Rees 
1971 lUgh tman & Shapi ro! [19771 ICo hn & Kul srudl 119781 
Syer & Ulmerll999tlMagorrian & Tremainel 19991) yields that 




the disruption rate continues to rise with decreasing energy for 
energies larger than the "critical energ)0", which in our case 



10000 



FIG . 19. — Disruption rate per star per log E, in units of Gyr , as a func- 
tion of energy. The rate drops quickly as the stellar distribution is depleted 
due to RR. 

is at E < 1. This is in accordance with Figure [19] As the 
figure shows, the rate drops quickly; this is due to the deple- 
tion of stars, and was als o found (though to a lesser ex tent, 
compare their Figure 4) in Hopman & Alexander (2006a). 

A limiting factor of our approach in this section is that the 
potential of our Galactic nucleus does not evolve; hence we 
cannot look at collective effects (e.g., instabilities). Ideally an 
iterative process should be used to self-consistently find the 
steady-state solution; we have shown that RR significantly 
depletes stars at inner radii. To what extent this will affect 
our steady-state solution is beyond the scope of the paper. We 
note the RR time does not formally depend on the stellar num- 
ber density until such energies and angular momenta that GR 
precession becomes important. 

7.2. A Depression in the Galactic Center 

Early studies of integrated starlight at the GC detected a dip 
in t he CO absorption strength within 15" of the MBH, Sgr 
A* dSellgren et al.| [T990; Hall er et al.ll 1996b . wh ich indicated 
a dec r ease in the density of old stars ( see also | Genzel et all 
[1991 12001 iFiger et all l2000L 12001 iSchodel et al.1 120071: 
Zhu et al.1 120081). Rec ent o bservations by iDo et alj (120091) . 
Buchholz etaiX ( 120091) . and iBartko et al J (120101) have con- 
firmed this suggestion and, in particular, revealed that the 
distribution of the late-type stars in the central region of the 
Galaxy is very different than expected for a relaxed density 
c usp aroun d an M BH. 

IDo et all d2009h use the number density profile of late-type 
giants to examine the structure of the nuclear star cluster 
in the innermost 0.16 pc of the GC. They find the surface 
stellar number density profile, oc i?~ r , is flat with 

T = 0.26 ± 0.24 and rule out all values of a > 1.0 at a 
confidence level of 99.7%. The slope measurement cannot 
constrain whether there is a hole in the stellar distribution 
or simply a very sha llow power law volume density profile. 
Buch holz et alj (120091) show that the late-type density func- 
tion is flat out to 10" with Y = -0.7 ± 0.09 and inverts in 
the inner 6", with power-law slope Y = 0.17 ± 0.09. They 
conclude that the stellar population in the innermost ~ 0.2 pc 
is depleted o f both bright and f aint giants. These results are 
confirmed bv lBartko et alj d2010l) who find that late-type stars 
with rrik < 15.5 have a flat distribution inside of 10". 



4 The critical energy is approximately the energy where the change in an- 



gular momentum per orbit equals the size of the loss cone. 
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This result has changed our perception of the nuclear star 
cluster at the center of our Galaxy. Stellar cusp forma- 
tion theory predicts a volume density profile n(r) oc r~ a 
of re laxed stars with a = 7/4 for a single mass popula- 
tion dBahcall & Wolfl[l976l) and has been confirmed in many 
theoretical p apers with different methods, including F okker- 
Planc k (e.gJCohn & Kulsrudll978tlMurphv et all 1991 . MC 
(e.g., IShapiro & Marchantlll978t IFreitag & BenzlhoOZ). and 
AT-body methods dPreto et all l2004t iBaumgardt et all l2004t 
iPreto & Amaro-Seoanel 12010). Theory also predicts that a 
multi-mass stellar cluster will mass segregate to a differential 
distribution with the more massive populations being more 
centrally concentrated. The indices of the power-law density 
profiles for the dif ferent mass populations can vary between 
3/2 < a < ll/4JBahcaU&Wolflfl977t IFreitag etal]l2006t 
I Alexander & Hopmanll2009l) . 

There have been two types of solutions put forward in the 
literature to explain the discrepancy between cusp formation 
theory and observations of the old stellar population in the 
GC. The first maintains that the stellar cusp in the GC is re- 
laxed as predicted, but additional physical mechanisms have 
depleted the old stars in the cent ral region. Such mechanisms 
include strong mass segregation dAlexander & Hopman 2009) 
and envelope de struction of giants by stellar collisions (see 
iDale et all 120091 and references within). The second solu- 
tion proposes an unrelaxed cusp, either by the depletion of 
the stars due to the infall of an intermediate- mass black hole 
(IMBH; Levin 2006: iBaumgardt et alj 120061) or from a bi - 
nary black hole merger dMerritt & Szellll2006tlMe"rritj|2010l) . 
Kinematical information for the old stellar population in the 
GC however shows a lack of eviden ce for any large-sca le dis- 
turbance of the o ld stellar cluster (Trippe et al. 2008]), and 



IShen et aT] (120101) find that the Galaxy shows no observational 
sign that it suffered a major merger in the past 9 — 10 Gyr. In 
addition to this, the dynamics of the S-stars and of the MBH 
greatly restricts the param eter space available for the p resence 
of an IMBH in the GC dGualandris & Merrittl 12009). For a 
summary see the latest review by Genzel et al.l ( 20 1 ) . 

Here we propose that these observations can be explained in 
part by a depression of stars carved out by RR acting together 
with tidal disruption. The underlying reason is that stars move 
rapidly in angular momentum space and enter the loss cone 
at a higher rate than are replaced by energy diffusion. We 
show that when we initialize the system as a cusp that, even 
though steady-state is reached in longer than a Hubble time, 
the depression forms on a shorter timescale. 

We perform simulations of stars, both for m = 1M Q and 
rn = 10M Q , initially distributed in a Bahcall-Wolf (BW) cusp 
around the MBH, such that dN/dE oc E~ 9 / 4 . We present 
first the results for m = 10Af Q , chosen to reflect the domi- 
nant stellar population at small radii. We plot the stellar dis- 
tribution function f{E), which scales as E x / A for BW cusp, 
at various times during the simulation as the stars move in an- 
gular momentum space due to RR; see Figure [20] In ~ 1 Gyr 
a depression is carved out of the population of stars, signif- 
icantly depleting the stellar population around the MBH out 
to ~ 0.2 pc. To conclusively demonstrate that RR is the rea- 
son for this depression, we also show the distribution of stars 
for simulations in which we have switched RR off (by set- 
ting 4>i = 0i = 0, o"i =^ 0) such that angular momentum 
relaxation follows a random walk. In this case no depression 
forms. These simulations show that, as the cusp is destroyed 
on timescales of a few Gyr due to RR, a BW cusp is not a so- 
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FIG. 20. — Carving out of a depression in the distribution of stars around 
an MBH for simulations in which m = 10-Mq. The broken orange line 
shows the stellar distribution at the start of the simulation (the solid black 
line above indicates a slope of 1/4). The thick colored lines show the results 
for RR while the corresponding (in color and style) thin lines show results for 
simulation in which RR is switched off. The energy, or equivalently radius, 
at which RR becomes efficient at carving out the depression is clearly seen 
on the scale of ~ 0. 1 pc. 




0.0001 



FIG. 21. — Number density of stars as a function of radius (in parsec and 
arcsecond) from the MBH for m = IOMq simulations. The stellar distri- 
bution begins to deviate from the BW solution at ~ 0.5 pc. The black lines 
indicate a slope of —1.75 (top) and —1.0 (bottom right). 



lution to the distribution of stars around an MBH. We see the 
same appearance of a depression when we substitute our non- 
resonant parameters (both energy and angular momentum) for 
those of EilonetalJ ( 120091) . 

In Figure |2T1 we plot the number density of stars as a func- 
tion of radius from the MBH. We do this by sampling each 
stellar orbit randomly in mean anomaly. Black lines indicate 
the initial BW number density slope of —1.75 and a slope 
of —1. The slope begins to deviate, albeit gently, from the 
BW solution at ~ 0.5 pc and decreases to —1 at ~ 0.1 pc. 
We present results from m — IMq simulations in Figures 
l22l and [23] even though, due to mass segregation, we do not 
expect low-mass stars to dominate the potential so close to 
an MBH. The major difference between the two simulations 
is the timescale on which the depression forms, 1 Gyr and 
10 Gyr, respectively. 

To reproduce the radial extent of the depletion of old 



14 



MADIGAN, HOPMAN AND LEVIN 



-'[pel 

0.01 0.001 



10 



a o.i 



0.01 



0.001 



t = 




1 Gyr 




5 Civ: 










A : A 



100 
log(E) 



1000 



10000 



FIG. 22. — Carving out of a depression in the distribution of stars around 
an MBH for simulations in which m = IMq. The timescale on which the 
depression develops is longer than for simulations in which the potential is 
built up of IOMq stars. Stars which have an energy E > 1000 survive 
longer than those with 100 > E > 1000 as RR is ineffective so close to the 
MBH due to GR. 
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FIG. 23. — Number density of stars as a function of radius (in parsec and 
arcsecond) from the MBH for m = IMq simulations. The depression in 
stars takes an order of magnitude longer in time to develop than the m = 
IOMq simulations. 



stars i n the observations of | Dp et al.l (120091) . iBuchholz et al.1 
(2009), and Ba rtko et al.l (120101) . the size of the hole in energy 
space must be significantly larger, by a factor of 5 — 10, 
than what we find in our simulations. In addition, a hole in 
energy space alone cannot give rise to an inverted profile in 
surface number density. We stress however that the potential 
in our model does not evolve in response to RR. As the 
depression develops and the number of stars at small radii 
decreases, the energy relaxation timescale, £nr oc N -1 , will 
increase. As N is not important for the RR timescale for all 
energies but the largest (GR precession at small radii affects 
the precession rate), this could further enlarge the depression 
as the rate of stars flowing inward toward the MBH decreases. 

In Figure [24] we demonstrate the increase in the NR 
timescale in response to the hole in stars. We use a fit to the 
stellar density profile at 1 Gyr from Figure [21] to recalculate 
£nr, normalizing the profile such that the mass in stars at the 
radius of influence is equal to that of the MBH. We plot the 



100000 



1000 




r [pc] 

FIG. 24. — NR (energy) timescale changing in response to the stellar po- 
tential (m = IOMq) as it develops from a BW cusp to one with a hole in 
energy space (from t^J^ p to t{^ c ). The RR timescale is plotted here without 
GR precession as M, /mP and, at large radii, does not depend on the number 
of stars. Hence this timescale should remain constant above ~ 10~ 2 pc as 
the NR timescale increases. The radius at which is equal to <jmr moves 
outwards in response to the changes and hence the outer radius of the hole in 
stars will increase. 



RR timescale without GR precession as relativistic effects are 
unimportant at radii greater than ~ 10~ 2 ; see Figure [15] As 
the outer extent of the hole at ~ 10 _1 pc is at a much larger 
distance from the MBH, the RR timescale at this scale will 
not be affected by changes in N. The increase in the radius at 
which <nr = £rr strongly suggests that the outer radius of 
the hole should increase in response to the change in stellar 
potential. 

One implication of this depression is that we expect 
no gravitational wav e bursts to be observed from the GC 
(Hop man et alj |2007). For the implications of the existence 
of a depletion of stars in the vicinity of an MBH on predicted 
event rates of gravitation al wave inspir als of compact objects 
(EMRIs) and bursts, see lMerrittl (120 lOt) . 

7.3. Dynamical evolution of the S-stars 

At a distance of ~ 0.01 pc from the MBH in the GC, there 
is a cluster of young stars known as the S-stars. This clus- 
ter has been the subject of extensive observation al and the- 
oretica l wo rk, for which we re fer the reader to lAlexanderl 
(2005) and iGenzel et alj ( 120101) . The stars are spectroscop- 
ically consistent with being late O- to early B-type dwarfs, 
implying that they have lifetimes tha t are limited to between 
< 6 and 400 Mvr jGhe7eTal]|2003t lEisenhauer et aT]|2005l: 
iMartins et alj 120081) . The upper constraint comes from the 
limit on the age of a B-type main-sequence star, while the 
lower constraint is a maximum age limit on a specific star 
(S2/S0-2). One particularly interesting question concerns the 
origin of these stars. There is a consensus in the literature 
in that they cannot have formed at their current location due 
to the tidal field of the MBH. Several scenarios incorporating 
the formation of these stars further out from the MBH, and 
transport in ward through dynamica l processes, have been pro- 
posed (e.g. [Gould & Ouilleril2003HAlexander & Livid | 2004[: 
lLevinl l2007l: iPerets et alJl2007HLockmann et alJl2008Tl2009l: 



Madigan et al. 20091: iMerritt et al.ll2009HGrivll2010l) . 

Many lead ing hypotheses involve binary disruptions dHillsl 
1988, 1991). In this scenario a close encounter between a 
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binary system and an MBH leads to an exchange interaction 
in which one of the stars is captured by the MBH, and the 
other becomes unbound and leaves the system with a high 
velocity. For capture, the periapsis r p of the binary's orbit 
(with combined mass mbm and semi-major axis abi n ) must 
come within the tidal radius 



/2M. 

\ m bin 



1/3 



Obi; 



(46) 



The captured star will remain with a semi-major axis a cap that 
scales as a cap ~ (M,/mbi n ) 2/ ' 3 <2bin7 and its eccentricity can 
be approximated as 1 — e ~ (mbm/M.) 1 / 3 . 

A unifying aspect of all proposed formation mechanisms is 
that none of them lead directly to orbital characteristics that 
are in agreement with those of the S-stars; in particular, the 
binary disruption mechanism leads to highly eccentric orbits. 
Since the local RR time is shorter than the maximum age of 
the S-stars, it has been suggested that following the arrival 
of the S-stars on tightly bound orbits to the MBH, a subse- 
quent phase of RR has driven them to their c urrent distribution 
(Hop man & Alexanderll2006atlLevinll2007l). RR evolution of 
S-stars has been studied by iPerets et al.l (120091) by means of 
direct A-body simulations. They find that high-eccentricity 
orbits can evolve within a short time (20 Myr), to an eccen- 
tricity distribution that is statistically indistinguishable to that 
of the S-stars. They conclude that the S-star orbital param- 
eters are consistent with a binary disruption formation sce- 
nario. However, these simulations do not include GR preces- 
sion, which we have found to be important for the angular 
momentum evolution of stars at these radii. 

7.3.1. Simulations 

We use our ARMA code, which includes GR precession, 
to study the evolution of the S-star orbital parameters within 
the binary disruption context. We initialize our test stars on 
high eccentricity orbits, e = 0.97 (m^n ~ 20M Q ), with 
initial semi-ma j or ax es corresponding to those reported by 
iGillessen et al.l (2009) for the S-stars. We assume a distance 
to SerA* of 8.0 kpc dEisenhauer et alJ2003HGhez et al.l2008b 
IGillessen et al .1120091) . For this distance, one arcsecond corre- 
sponds to a distance of 0.0388 pc. We do not include S-stars 
that are associated with the disk of OAVR stars at larger radii 
for two reasons. The first is that they most likely originated 
from the disk with lower initial orbital eccentricities than we 
simulate here. Second, RR is not the most effective mecha- 
nism for angular momentum evolution at their radii; torques 
from stars at the inner edge of the disk will dominate. We 
run the simulation one hundred times for each test star to look 
statistically at the resulting distribution. 

We simulate two different initial setups. The first, which 
we call the continuous scenario, assumes that the tidal disrup- 
tion of binaries is a n ongoing process, as is most likely in the 
case considered by IPerets et ail J200 7). where individual bi- 
naries are scattered into the loss cone by massive perturbers. 
The same initial conditions also hold if low angular momen- 
tum binaries are generated as a result of triaxiality of the sur- 
rounding potential. For this study we start evolving the stellar 
orbits at times that are randomly distributed between [0, t max ], 
and continue the simulation until the time t = t max . In our 
second setup, which we call the burst scenario, we consider 
a situation in which S-stars are formed at once on eccentric 
orbits at t = 0, and follow their evolution until t = t max \ this 
enables us to compare directly to the results of Pere ts et al.l 



(2009) in whose simulations all stars commence their evolu- 
tion simultaneously. A burst of S-star formation, over a short 
period of a few Myr, could for exa mple result from an insta- 
bility in an eccentric stellar disk (Madigan et al.l 120091) : we 
note t hat the instabili ty is even more effective for a depressed 
cusp (Madiga nl2010l) . 

We use two different models for the galactic nucleus. For 
each we re-compute the ARMA parameters (<f>, 9, cr) and 
tNR for the new profile and re-run the MC simulations using 
Equations d40l > and (HTt with the updated values. Model A 
is similar to the galactic nucleus model discussed in Section 
|2] The stellar cusp has a density profile n(r) oc r~ a with 
a = 1.75, and the number of stars is normalized with 
N(< 0.01 pc) = 222 and m = 10. We chose this model 
such that it faithfully reflects the dominant stellar content 
at a distance of 0.01 pc, where RR is most important and 
where the S-stars are located. We use Model B to simulate 
the response of the S-stars to a depression in the GC, as 
found in the previous section. This model has a smaller 
density power-law index of a = 0.5 (which is within the 
90% bound aries of t he bo otstrapped values for a in the GC 
as found by iMerritj (120101) ') and N(< 0.01 pc) = 148. Less 
mass at small radii leads to slower mass precession, pushing 
to a lower eccentricity the value at which mass precession 
cancels with GR precession; see Figure [l6]for an illustration 
of this effect where a = 1.75. Interestingly for all values of 
N(< 0.01 pc) < 120, this eccentricity lies at ~ 0.7 a location 
in eccentricity space where we do not observe S-stars, which 
causes a flattening of the observed cumulative distribution. 



TABLE 3 
S-star Parameters 



Model A 


imax (Myr) 


6 


10 


20 


100 


p a 


Burst d 


0.012 


0.147 


0.007 


5e-4 




Continuous 6 


2e-4 


0.002 


0.05 


0.02 


f > 0.97 6 


Burst 


.40 


.29 


.24 


.07 




Continuous 


.47 


.40 


.33 


.16 


DR C 


Burst 


.23 


.32 


.46 


.72 




Continuous 


.09 


.17 


.29 


.55 


Model B 


p a 


Burst 


0.002 


0.04 


0.100 


0.0005 




Continuous 


1.5e-5 


0.0002 


0.006 


0.06 


f > 0.97 b 


Burst 


.45 


.35 


.22 


.07 




Continuous 


.50 


.47 


.38 


.19 


DR C 


Burst 


.17 


.29 


.45 


.70 




Continuous 


.08 


.16 


.26 


.55 



a Two-dimensional (a, e) Kolmogorov-Smirnov p-values. 
b Fraction of S-stars with e > 0.97 at end of simulation. 



c Disruption rate (fraction) of S-stars. 

d Stars initialized at t = 0. 

e Stars randomly initialized between [0, 



7.3.2. Results 



For both models we find the cumulative eccentricity distri- 
bution of the stars at diffe rent times and compar e them to the 
observed distribution from Gill essen et al.l(| 2009): see Figures 
[25] and [26] In Table [3] we present results from nonparametric 
two-dimensional Kolmogorov-Smirnov (2DKS) testing be- 
tween the observed cumulative distribution and the simula- 
tions. The two-sample KS test checks whether two data sam- 
ples come from different distributions (low p values). We also 
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give the fraction of stars that have an eccentricity e > 0.97 at 
the end of the simulation. The total disruption rate (DR) gives 
an indication of how many progenitors would be needed to 
make up the population of S-stars observed today. 

Model A: The best fit to the S-star observations are for the 
burst scenario with i max = 6 — lOMyr. Greater values of 
imax produce too many low eccentricity stars in our simula- 
tions; as the RR time is longer at these eccentricities stars will 
tend to remain at these low values. In the continuous scenario 
the simulations conclude with many more high eccentricity 
stars as they do not have as much time to move away from 
their original eccentricities. Hence the best-fit result for this 
scenario is i max = 10 — 20 Myr. In both cases the 2DKS p- 
values <C 1, suggesting that no simulated distribution matches 
very well the observed. Both scenarios predict a large number 
(7 — 45%) of high eccentricity stars with e > 0.97. 

Model B: The best match to the observations is the burst 
scenario with t max = 20 Myr, though the continuous scenario 
has a relatively high p- value at a time of £ max = 100 Myr. In 
both cases a longer t max is preferred as this model has less 
mass, and hence GR precession will dominate the precession 
rate of these stars down to lower eccentricities, rendering RR 
less effective at high eccentricities. Again the p-values are 
low in most cases due to the excess of high eccentricity orbits 
relative to the observations. 

For both models, simulations with i max > 100 Myr pro- 
duce too many low eccentricity stars with respect to the ob- 
servations. Simulations with t max = 200 — 400 Myr have in- 
creasingly smaller p- values for this reason. The small number 
of observed low eccentricity S-stars could in principle con- 
strain the time S-stars have spent at these radii (their "dynam- 
ical lifetimes"). Although the S-stars can in theory have ages 
up to ~ 400 Myr, we find that, in the binary disruption con- 
text, their dynamical lifet imes at thes e radii are < 100 Myr. 
Furthermore, as noted by Per ets et al.l ((2009), the lack of ob- 
served low eccentricity S-stars is difficult to reconcile with 
formation scenarios which bring the stars to these radii on 
initially low eccentricity orbits. The RR timescale t„ is long 
and stars will retain their initial values for long times. 

For all values of i max there is a clear discrepancy between 
the lack of high eccentricity S-stars orbits observed and the 
results from both formation scenarios and models. To fur- 
ther illustrate this problem we show a scatter plot from our 
simulations in Figure [27] of eccentricity e as a function of 
semi-major axis a after 6 Myr. We over-plot the observed 
S-star parameters, indicating separately those associated with 
the disk. It is clear that while the simulations can reproduce 
the range of (a, e) in which the S-stars are observed, our the- 
ory over-predicts the amount of high eccentricity orbits with 
respect to the observations. Such orbits at these radii are ex- 
pected; stars with high orbital eccentricities, although they 
experience large torques (Equation (l30j>), have short coher- 
ence times (Equation ( |29l ) and hence long RR times (Equa- 
tion (|3TT>). Consequently, their angular momentum evolution 
is sluggish. At any given time, we would expect to see a 
population of stars in this region. Therefore, we find it dif- 
ficult to explain why so few high eccentricity S-stars have 
been observed; the m aximum is e = 0.963 ± 0.006 (S14; 
Gillessen et al1 l2009): e = 0.974 ± 0.016 (S0-16: lGhezetaD 
20051) . 

An interesting solution to this problem is that there are S- 
stars with higher orbital eccentricities but that their orbital el- 
ements have not yet been determined due to the observational 
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FIG. 25. — Cumulative distribution of eccentricities for the burst scenario 
in which the stars' initial time is zero. The red step function is the ob- 
se rved cumulative eccen tricity distribution of S-stars in the GC as found 
bv lGillessen"et~aTl f2009h . The dashed orange line is a thermal distribution 
N(< e) = e 2 . The other lines are results from MC simulations, where the 
evolution time t max and eccentricity e are indicated in the legend. Thick 
(thin) lines are the results for Model A (B). 
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FIG. 26. — Cumulative distribution of eccentricities for the continuous 
scenario in which the stars' initial time t is randomly distributed between 
[0, tma X ]- The red step function i s the observed cumula tive eccentricity dis- 
tribution of S-stars in the GC from Gillessen et al. (2009). The dashed orange 
line is a thermal distribution N(< e) = e 2 . The other lines are results from 
MC simulations, where the maximum evolution time t max and initial ec- 
centricity e are indicated in the legend. Thick (thin) lines are the results for 
Model A (B). 
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FIG. 27. — Scatter plot of eccentricity, e, of the S-stars, as a function of 
semi-major axis, a, on a log scale for the Model B burst scenario. The small 
orange dots represent the results of our sim ulations after 6 Myr, the blue 
squares are data from Gillessen et al. (2009) and the black open circles are 
those S-stars consistent with being members of the CW disk. 



detection limits for stellar accelerations. ISchodel et alJ (120031) 
show that there is a bias against detecting high eccentricity 
orbits at the average S-star radius of 0.01 pc (see their Fig- 
ure 14). Future observations will be able to reveal if this can 
account for the discrepancy (Weinberg et al. 2005). The so- 
lution on the other hand may involve a new stellar dynam- 
ical me chanism which is n ot accounted for in our ARMA 
model. iMerritt et al.1 (1201 ll) have recently discovered a new 
dynamical barrier in angular momentum space, termed the 
"Schwarzschild barrier", which can reflect stars orbiting close 
to the MBH from very high eccentricity orbits (e > 0.99) to 
less eccentric ones. This could well explain why not so many 
S-stars have been observed on highly eccentric orbits as RR 
alone predicts. 

8. SUMMARY 

We have used several techniques to study the angular mo- 
mentum evolution of stars orbiting MBHs, with a focus on 
secular effects. First, we carried out A-body simulations to 
generate time series of angular momentum changes for stars 
with various initial eccentricities. We quantified the evolution 
in an ARMA( 1,1) model. We then used this model to gener- 
ate new, much-longer time series in our MC code, which en- 
abled us to study the steady-state distribution of stellar orbits 
around an MBH, and the evolution of the angular momentum 
distribution of possible progenitors of the S-stars. 

We have shown that the steady-state angular momentum 
distribution of stars around an MBH is not isotropic when RR 
is a dominant relaxation mechanism. For the GC this corre- 
spo nds to the innermost couple of tenths of parsecs. We note 
that lSchodel et al.l (|2009) find a slight degree of anisotropy in 
the velocity distribution of late-type stars in the GC within 
the innermost 6" (~ 0.2 pc). It is of particular interest to 



note that close to the MBH, we find that the angular momen- 
tum distribution is steeper than linear, w hich may have conse- 
quences for the stability o f the system. iTremainel (120051) and 
Polyachen ko et al.1 (12007) showed that if (1) the distribution 
rises faster than linearly, and (2) the distribution of angular 
momenta vanishes at J = 0, the system is only neutrally sta- 
ble (or even unstable if it is flattened). The second condition is 
also satisfied, because of the presence of the loss cone: N = 
for J < Ji c . In a later paper, Polyachen ko et al.l (|2008) found 
that a non-monotonic distribution of angular momenta is re- 
quired for the system to develop a "gravitational loss cone 
instability". This condition is also satisfied for a range of en- 
ergies. An assessment of the importance and consequences of 
these instabilities is beyond the scope of this work. 

We have furthermore demonstrated the presence of a de- 
pression in the distribution of stars near MBHs due to resonant 
tidal disruption. This depression forms on a timescale shorter 
than the Hubble time, even though a steady-state solution may 
not yet be reached. While it is tempting to invoke this mech- 
anism to explain the "hole" in late-type stars in the GC, we 
stress that our fiducial galactic nucleus model, a single-mass 
BW distribution of stars around an MBH, is greatly simplified 
compared to reality and produces an outer hole radius which 
is too small to explain the observations. We expect the GC to 
be far more complex, with ongoing star formation, multi-mass 
stellar populations, mass segregation, physical collisions be- 
tween stars and additional sources of relaxation (e.g., the cir- 
cumnuclearring, the CW disk). Nevertheless, there is merit to 
our approach: even with our simplified picture, we illustrate 
the importance of RR on stellar dynamics close to an MBH. 

Finally we have shown that the (a, e) distribution of the S- 
stars does not match that expected by the binary disruption 
scenario in which the stars start their lives on highly eccen- 
tric orbits and evolve under the RR mechanism. We surmise 
that there are a number of high eccentricity S-stars whose 
orbital parameters have not yet been derived , though the 
"Schw arzschild barrier" recently discovered by Merr itt et al.l 
may provided a dynami cal solution t o this problem. 
We have confirmed the result bv lPerets et al.l (120091) in that it 
is unlikely that the S-stars were formed on low eccentricity 
orbits. We place a constraint on the dynamical age of the S- 
stars of 100 Myr from the small number of low eccentricity 
orbits observed. 
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APPENDIX 
A. DESCRIPTION OF iV-BODY CODE 



We have developed a new integrator, specifically designed for stellar dynamical simulations in near- Keplerian potentials. 
A dis tinctive feature of this code is the separation between "test" particles and "field" particles (see, e.g., iRauch & Tremainel 
1996). The field particles are n ot true A-body particles; they move on Keplerian orbits which precess according to an analytic 
prescription (see Equation (IC8I 1). and do not directly react to the gravitational potential from their neighbors. Conversely, the test 
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FIG. 28. — Fractional difference between the theoretical value for apsidal precession and that achieved by our integrator for a test star of eccentricity e = 0.99 
as a function of time in units of the precession time of the test-star. The angle through which the test star has precessed is given by <j>. The analytical formula is 
adjusted for the test star's wandering in semi-major axis. 



particles are true A-body particles, moving on Keplerian orbits and precessing due to the gravitational potential of all the other 
particles in the system. Both test and field particles' orbits precess due to GR effects. This reduces the force calculation from 
0{N 2 ) — > 0(nN), where A is the number of field particles and n is the number of test particles in the simulation. 

The main difference between our code and many others, is that instead of considering perturbations on a star moving on a 
straight line with constant velocity, we consider perturbations on a star moving on a precessing Kepler ellipse. Our algori thm 
is based on a mixed-variable symplectic method ( Wisdom & Hormani ri99U [Kinoshita et alJll99Tt ISaha & Tremaindl 19921) . so 
called due to the switching of the coordinate system from Cartesian (in which the perturbations to the stellar orbit are calculated) 
to one based on Kepler elements (to determine the influence of central Keplerian force). These forms of symplectic integrator 
split the Hamiltonian, H, of the system into individually integrable parts, namely a dominant Keplerian part due to the central 
object, i.e., MBH, and a smaller perturbation part from the surrounding particles, i.e., stellar cluster, 



The algorithm incorporates the well-known leapfrog scheme where a single time step (r) is composed of three components, 
kick (r/2)-drift (r)-kick (r/2). The kick stage corresponds to ^interaction and produces a change of momenta in fixed Cartesian 
coordinates. The drift stage describes motion under i?Kcpior and produces movement along unperturbed Kepler ellipses. To 
integrate along these ellipses, we exploit the heavily-studied solutions to Kepler's equation. This is a less computationally efficient 
method than direct integration, but accurate to machine precision. This is of utmost importance in simulations involving RR as 
spurious precession due to build-up in orbital integration errors will lead to erroneous results. The computation is performed in 
Kepler elements, using Gauss' / and g functions. We use universal variables which allow the particle to move sm oothly betwee n 
elliptical, hyperbolic, and parabolic orbits, making use of Stumpff function s. For details w e refer the reader to iDanbyl (f 1992). 
For the simulations presented in this paper, we use a fourth-order integrator ( Yoshida 1990) which is achieved by concatenating 
second order leapfrog steps in the ratio (2 - 2 1 / 3 )" 1 : -2 1 / 3 (2 - 2 1 / 3 )" 1 : (2 - 2 1 / 3 )" 1 . 

We use adaptive time stepping to resolve high eccentricity orbits at periapsis, and to accurately treat close encounters between 
particles, hence our integrator is not symplectic. To reduce the energy errors incurred with the loss of symplecticity, we 
implement time symmetry in the algorithm. These steps must be sufficient to avoid artificial apsidal p recession of stellar orbits , 
particularly at high eccentricity where Wisdom-Holman integrators are known to perform poorly ( Rauch & Holman 1999). 
To demonstrate this we compare the rate of change, as a function of time, of the direction of the eccentricity vector of a test 
star using our A-body integrator with our analytical prediction. Figure [28] shows our results in terms of the angle through 
whi ch the test star has precessed, <f>, for e = 0.99, a = 0.01 pc orbits in the gravitational potential as described in Section 
12.11 We take the average of nine simulations for half a precession time of the test star. Extrapolating to a full precession 
time we find a fractional difference of S(p/<p = 0.00148 between our analytical estimate and that of the test star precession 
rate. This verifies that the long RR timescale at high eccentricities, as shown in Figure [3] is a result of attested rapid apsidal 
precession. To further demonstrate the accuracy of the integrator with respect to the apsidal precession rate, we run additional 
calibration simulations for high-eccentricity orbits (e = 0.96, 0.98) in order to compare their angular momentum relaxation 
timescale (calculated using Equation ([5])) with the formula we derive for RR in Equation PIT ). Figure [29] demonstrates the 
agreement between increase of the RR timescale with eccentricity as retrieved from the simulations and that expected from theory. 

An additional feature of the code is the option to implement "reflective" walls, which bounce field particles back into the sphere 
of integration. Without these walls, particles on eccentric orbits nearing apoapsis move outside the sphere of interest and alter the 
density profile of the background cluster of particles, and hence the (stellar) potential. With the walls "on", a particle reaching the 
reflective boundary is reflected directly back but with a different random initial location on the sphere, to avoid generating fixed 
orbital arcs in the system which would artificially enhance RR. The reflective walls are most useful in situations where we want 
to examine how a test particle reacts in a defined environment with a specific background (stellar) potential. The prescription 
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FIG. 29. — Timescale for angular momentum J changes of the order of the circular angular momentum for massless test stars in our TV-body simulations, 
presented on a log scale of time in units of Myr; see Equation (6). We tak e se veral At values to get order-of-magnitude estimates for this timescale. Over-plotted 
in red is the theoretical expectation for this timescale based on Equation )3U . 



we use for the time step (r) is based on orbital phase, the collisional timescale between p articles and fre e-fall timescale between 
particles. We treat gravitational softening between particles with the compact K2 kernel dDehnenl2001l) . 
We have parallelized our code for shared-memory architecture on single multiple-core processors using OPENMR 



B. ENERGY EVOLUTION AND CUSP FORMATION 



Many different definitions of relaxation times are used in the literature. In Section 12.21 we consider the energy diffusion 
timescale ts, defined as ts = E 2 / Dee where Dee = ((AE) 2 ) / At is the energy diffusion coefficient. This is the timescale 
for order unity steps in energy. In order to be explicit, we now compare our timescale (see Figure [3] and Equation (|5]l) with a 
commonly used reference time for relaxation processes defined in Spitzej] dl987l) : 



■M 



3((A W ||)2)^ A 



(Bl) 



= 0.065- 



M 



(Gm) 2 n In A' 

We note that this is equivalent to the timescale given by iBinney & Tremaind ((2008) for a Maxwellian velocity distribution, 
v m = VScr, where Vm{i") is me mean velocity of a star and <j(r) is the one-dimensional velocity dispersion of the stellar 
distribution at radius r, 



0.34cr J 



(Gm) 2 n In A' 



(B2) 



The velocity dispersion within a stellar cusp with power-law density profile n(r) oc r a can be calculated from the Jeans equation 
for spherical isotropic systems as 



a 2 (r) ee - 
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(B3) 



where M(r) is the mass within a radius r and v c is the circular velocity. For the comparison of t E and t r we take Keplerian 
parameters for the circular velocity at a given radius such that 



o 2 = 



1 GM. 



(B4) 



( I Alexa nder 2003), and set the Coulomb logari thm to A = M. /m (Ba hcall & Wolfll976l) . For the cusp model used in our A-body 
simulations we find from Equations d23l and iB2\ 



E 2 (M. 

t E = jt— = 0.26 — 
Dee \ m 



1 1 
A< In A 



P = 2.18 t r 



(B5) 



with tf = 156 Myr, ts = 340 Myr for r = 0.01 pc. Thus with this definition, the energy of a star changes by order unity in 
roughly twice a relaxation time. 



We compare our result for the energy relaxation time to that found by lEilon et ail (2009) using an entirely different inte- 
grator. Their A-body simulat ions incorporate a fifth-order Runge-Kutta algorithm with optional pairwise KS regularization 
dKusta anheimo & Stiefell 119651) , without the need for gravitational softening. We use their numerically derived value (see 
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FIG. 30. — Stellar distribution function f(E) for single mass stars as a function of energy E for different times t = (1 — 10) tg, where is the Keplerian 
energy diffusion timescale evaluated at E = 1. steady-state is reached within 10 tE, the solution of which is a BW cusp with a = 7/4, represented by the 
dotted green line. 



their Table 1) of (a) = (aA/Vln A), averaged over all stars in the cluster, and their expression for the energy relaxation time 
Tjfpj = (M,/m) 2 P/(NajJ, and compare with our Equation (f23t with its own numerically derived parameter. We find an 
energy relaxation timescale which is a factor of ~ four longer than in their paper (i.e., their energy diffusion timescale is about 
half of t r ). Our maximum radius at 0.03 pc means we do not simulate several "octaves" in energy space which could contribute 
to energy relaxation. In a Keplerian cusp, the variance in the velocity changes of a star, (Av) 2 , is not independent of radius r but 
scales with the power-law density index a. In our simulations a — 7/4, and hence the number of octaves from the cutoff radius, 
r = 0.03 pc, out to the radius of influence, rh — 2.31 pc, can account for at most a factor of ~ two difference in tE- The inner 
cutoff radius at r = 0.003 pc contributes a smaller factor of difference to the energy diffusion timescale. 

We now consider the timescale for building a cusp due to energy diffusion. For this purpose, we use our MC code without 
angular momentum evolution. In Figure|30] we show the distribution function for several times, where the system was initialized 
such that all stars start at E = 1, corresponding to half the radius of influence. From this figure, we find that it takes about 
10 tE(E = 1, r = 0.5 r/i ) or 8 ^^(r = r/,) to form a cusp. The result that it takes much longer than tE is also implicitly found 
in iBahcall & Wohl (119761) . They give an expression for Dee (called C2(E, t) in their paper); evaluating the pre-factor in their 
expression due to only the contribution of bou nd stars in steady-state, one finds that c^iE = 1) = 16/t^^ p , where is the 
reference time in which Bahc all & Wolf) (119761) find a cusp is formed. 

iBaumgardt et al.l (|200 4) a nd iPreto et alj d2004) were the first to show the formatio n of a cusp in dire ct A-body simulations. 
We use Equations (1B21 I and dB51 l to calculate the cusp formation times in Table 1 of IPreto et al.l (120041) in terms of our energy 
diffusion time. We find that in their simulations a cusp forms on a times cale of the order o f 10 tE- We conclude that energy 
diffusion proceeds in their simulations at a rate similar to ours. However. IPreto et all d2004) state that a cusp forms within one 
relaxation time. These two statements, apparently contradictory, are reconciled in two ways. First. IPreto et al.l (120041) define the 
relaxation time at a radius such that the mass in stars is equal to twice that of the MBH. Second, they include the non-Keplerian 
potential in their calculation of the relaxation time. It is valid to assume Keplerian parameters for the stellar velocity dispersion 
at the distances we are interested in close to the MBH. However, neglecting the contribution of the surrounding stellar mass leads 
to a significant numerical difference in the calculation of the relaxation time near the radius of influence compared with those 
which incorporate the non-Keplerian potential. 

To verify this, and to reconcile our conclusions, we adapt our definition of the energy diffusion time by incorporating the 
non-Keplerian potential, 



where A = 2 a 3 Mf l /M„ is the mass inside the radius of influence r^, t' E is the new energy diffusion timescale and we have 
used the relation E = rh/2a. The new energy diffusion coefficient is 



In Figure [3J] we show the results of MC simulations in which the new energy diffusion coefficient is used. The system was 
again initialized such that all stars start at E = 1. A BW cusp (a = 7/4) forms in about one relaxation time evaluated at 
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FIG. 31. — Stellar distribution function f(E) for single mass stars as a function of energy E for different times t = (0.3 — 1) t' E , where t'„ is the energy 
diffusion timescale which incorporates the non-Keplerian stellar potential and is evaluated at E = 0.287, r = 4.023 pc. The stars start in a delta function at 
E = 1 (and hence generate an artificial peak at this position) and form a cusp within one relaxation time 1 t' E . 



E = 0.287, the energy at which M/, = 2M # . This result emphasizes the importance of taking the stellar potential into account 
for the relaxation timescale. 

For a simplistic single mass (m = 1 M©) GC model (see Section 12. U . we find that at the radius of influence the Keplerian 
energy diffusion timescale is = 12 Gyr and a cusp forms in ~ 100 Gyr. Alternatively we find from Equation dB7l ) a cusp 
formation time of t' E ~ 90 Gyr. We stress that the time for cusp formation is much longer than the Hubble time, even if we 
halve it to account for the outer cutoff in radius in our TV-body simulations which suppresses the energy diffusion times cale. 
This is in contradiction to w hat is often claimed; one cause of the discrepancy is that in, for example, Bahcall & Wolf ( 1976) and 
Light man & Shapirol (119771) . energy diffusion proceeds at a rate exceeding the relaxation rate by a factor > 10, in contrast to what 
we find in our iV-body simulations. We note here that if energy diffusion did in fact proceed at this higher rate, then the energy 
diffusion timescale for the S-stars in the GC would be on the order of ~ 20 Myr which would have prof ound conseq uences for the 
(E, J) evolution of the cluster. It thus appears that the GC cannot have reach ed a steady-state (see a lso Merritt 2010), e ven though 
mass- segregation enhances the rate of cusp formation by a factor of a few dPreto & Amaro-Seoandl2010t iHopman & Madigan 
l2010h . 

C. PRECESSION DUE TO POWER-LAW STELLAR CUSP 

The presence of a nuclear stellar cluster in the vicinity of an MBH adds a small correction 5U(r) to the Keplerian gravitational 
potential energy U = —a/r of the system. As a consequence, paths of finite motion (i.e., orbits) are no longer closed, and with 
each revolution the periapsis of a star is displaced through a small angle 8(f). Given a power-law density profile of the stellar 
cluster n(r) oc r~ a , we can calculate this angle using the formula 

5<I >= JL(1 f r 2 5Ud(j)] , (CI) 



dL \L . 

where L = \J GM m a (1 — e 2 ) is the unnormalized angular momentum of a Keplerian orbit (lLandau & Lifshitz1ll969l) . We 
continue to use specific parameters in the following equations. The stellar mass within radius r is 

\ 3-a 

r 



M(r) = Mo l^-j , (C2) 

Mo, ro being proportionality constants which we normalize to the radius of influence, ro = r^, Ma = M.. Thus the potential 
felt by a massless particle due to this enclosed mass is 



5U = - / F(r')dr' 
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Inserting SU into Equation (ICU yields 
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We apply a change of coordinates here, replacing L with normalized angular momentum J = \/l — e 2 , such that L = ^/GM,aJ. 
Expressing radius r in terms of Kepler elements, semi-major axis a, eccentricity e and angle from position of periapsis (j>, 



r = a- = a- 

1 + e cos <p 



J 



l + x/l^J 2 . 



(C5) 



yields 



(2 -a) \r h 
2 



3-a 



(2 -a) 



f(e,a) 



f(e,a) 
M. 



(C6) 



where m is the mass of a single star, N < is the number enclosed within radius r and 
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The function /(e, a) is calculated numerically, and fit, for a given a. This returns a cluster precession time of 
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where P(a) is the period of an orbit with semi-major axis a. 

D. EQUATIONS OF ARMA( 1,1) MODEL 

We will use Equations ([H} and (O in our calculations; for clarity we repeat them here, dropping the label "1", 

AJ t = 0AJ t _i +0e t _i +e t , 
(e) = 0; (e t e s ) = a 2 8 t , s . 
D.l. Variance 



(Dl) 
(D2) 



(AJ t 2 ) = 2 (AJ t 2 _!) + 200<AJ t _iet_i) + 20(AJ t _ l£t ) 
+ tf 2 ( e 2 _ 1 }+20(e t _ 1 e t ) + ( e 2 > 

Using (A J 2 ) = (AJ^j), expanding A J t _i in the second and third terms, and applying (ete s ) = <T 2 5 t . s yields 

(AJ 2 )(l-0 2 ) = a 2 (2<b8 + 6 2 + l), (D4) 

which returns Equation ( flOl i: 

/A 9l l + 6» 2 + 26W> „ 

(AJt } = l-p a - (D5) 

D.2. Autocorrelation function 
The ACF for the ARMA model as described by Equation (ID 11 1, is defined as 

( AJ s+t A J : 
Expanding the numerator gives 

(AJ s+t AJ s ) = ^{AJ, +t _iAJ,_i) + 00<AJ, +t _ie,_;i) 
+ <£(A J s+t -ie s ) + ^<AJ,_ie H . t _i) 
+ </)(AJ s _ie s+t ) + 6» 2 (e s+t _ie s _i) 
+ 9(e s+t -ie s ) + 6(e s+t e s -i) + (e s e s+t ), 

where (AJ s+t AJ s ) = (j) 2 (A J s +t-i A J s _i). Recursively expanding the A J terms, and again using (e t € s ) — a 2 5t, s , reduces the 
fourth and fifth terms to zero, and greatly simplifies the expression to 

(aj , +iaj ,) . + (D8) 



Pi = ' T" (t>0). <D6) 
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Normalizing by the variance returns Equation ( fTTT i: 
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D.3. Variance at coherence time t$ 
To calculate the variance at the coherence time, (A Jj), we begin by calculating the variance after some time t, 

ft ft 
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(D9) 



{ AJ?) I [ P( tl -t 2 )dtidt 2 . 



(D10) 



Here pc^-t^) is the ACF of the torque at {t\ — £2) > 0, normalized with (A Jf) to be comparable with Equation ( fTTT i, 



P(U-t 2 ) 



From Equation (1321 we can write 



h(ti-t 2 ) 



(ti-ta) 



exp 



1 + (0 + 0) 2 /(l -^ 2 ) 
(*l-*2) 



Inserting this into Equation (IDIOK writing 



A = (AJ, 2 ) 



and taking an upper limit of t — t$ yields 



(AJ2)=A / dtx I exp 



1 + (0 + 0) 2 /(l ~<t> 2 ) 
(*1 - t 2 ) 



(Dll) 



(D12) 



(D13) 



o jo 



dto 



(e + ct>)(d + i/<f>) 



p 



(D14) 



where we have approximated (e + 1/e — 2) « 1. 
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